/*===================================================================
 * QuadElement.hpp
 * This is the 2D NURBS element data structure.
 * We define a class of QuadElement. It should have the following 
 * properties:
 * (1) The local basis function number and the mapped global numbering;
 * (2) Quadrature points and weights;
 * (3) Control points --> weights array (for NURBS construction);
 * (4) Element index;
 * (5) Evaluate the basis function and its derivatives up to 2nd order;
 * (6) Jacobian;
 * 
 * ------------------------------------------------------------------
 *  Author: Ju Liu, Sept 22nd 2011
 *          The University of Texas at Austin
 * ================================================================*/

#ifndef QUADELEMENT_HPP
#define QUADELEMENT_HPP

#include <vector>
#include <NURBS_Bezier.hpp>
#include <MeshInfoProcess2D.hpp>

class QuadElement
{
  public:
    QuadElement( int element_index,
        const class MeshInfo_2d &mesh2d,
        std::vector<double> ctrlPts,
        int cpdim,
        std::vector<NURBS_T::BezierElem*> seg_x,
        std::vector<NURBS_T::BezierElem*> seg_y );
    
    
    ~QuadElement();
    
    int enindex(){return en_index;}
    
    /* Evaluate the value of basis function and its 
     * derivative on quadrature points qua_x, qua_y.
     * The output should give R, dR_dx, dR_dy d2R_dx2
     * d2R_dy2 d2R_dxdy and Jacobi                        */
    void ElemShapeRoutine(double qua_x, double qua_y);
    
    double * R;
    double * dR_dx;
    double * dR_dy;
    double * d2R_dx;
    double * d2R_dy;
    double * d2R_dxdy;
    double Jacobi;


  private:
    /* Element basic information */
    int en_index;       // element global index
    int en_index_x;     // element index in x coordinate
    int en_index_y;     // element index in y coordinate
    
    int enLocBas;
    
    double en_hx;       // element hx
    double en_hy;       // element hy

    int * localIEN;     // local IEN array
    int sDegree, tDegree;            // order in x y direction
    int enqp_x, enqp_y; // number of element quadrature points
    
    double * qp_x;      // qua pts in x dir
    double * wq_x;      // qua weights in x dir
    double * qp_y;      // qua pts in y dir
    double * wq_y;      // qua weights in y dir
    
    double * ctrl_x;    // control points x coor
    double * ctrl_y;    // control points y coor
    double * ctrl_w;    // weights 

    double ** c_x;     // extraction operator in x (xi) direction
    double ** c_y;     // extraction operator in y (eta) direction

    void LinearSolver(double a, double b, double c, double d,
                     double f1, double f2, double f3, 
                     double &sol1, double &sol2, double &sol3);

};

#endif
